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Abstract 

The atmospheric muon flux is generated with CORSIKA 6.030 (using QGSJET high-energy interac- 
tion model) and fitted with the formula from Gaisser (1990). Spectra and mass distribution of primaries 
were taken according to the poli-gonato model of Horandel (2003). A convenient parameterization of 
the cos{9) cos{6*) correction, originally tabulated in Volkova (1969), is presented. This correction, 
together with muon energy losses and decay, is shown to significantly improve the fit to the simulated data. 

o : 1 Introduction 

^ The zenith-angle dependent energy distribution of muons, generated in showers initiated by high-energy 
cosmic rays in the atmosphere, can be conveniently parametrized [T| as 

dN _ 0.14 /^_^V^ ( ^ 0.054 \ 

l> . ~ cm2 sr s GeV ' ' 1 GeV / ' 1 i + i-i^a^cq^W + i , i-I-^mcq^^W I ^ ^ 

^ ^ \ ^ 115 GeV ^ 850 GcV / 

^ This formula is often used in the data analyses of the underground muon and neutrino detectors, e.g., for 
00 the vertical depth-intensity relation calculation or the estimates of the charm component contribution to the 
^ muon flux. The relative normalization A and the spectral index 7 were varied in this work to fit the results of 
] CORSIKA [2 1 simulations. With the assumption of scaling in the high-energy hadron-nucleus interactions 
Q ' used in derivation of this formula, the muon flux spectral index 7 should be the same as that for the 
■ primaries. However, in the parameterization of the cosmic ray energy spectrum and mass distribution used 
^ here [4J, different primaries have different spectral indices. Therefore the expected value of the muon flux 
Q-i spectral index 7 is not immediately obvious, but should be close to some "weighted average" of spectral 
indices of the individual cosmic ray components. In |5| values of tt and K critical energies |6| at = and 
^ the ratio of K to vr-generated muons were also varied. In this work varying these parameters did not improve 
the fits, so they have been fixed at 115 GeV, 850 GeV, and 0.054 as indicated in the formula above (as in I.1J). 

According to B, the cosmic ray energy spectrum can be described by a sum of components with Z=l-92, 
each given by 

E 

Ez , 

The best description of cosmic ray data within this model is achieved using rigidity-dependent cutoff energies 
Ez = Ep-Z. Values of Ep, A7, and Cc, as well as a table of values of <l>^ and 7^ used in this work are given in 
||4|. Since CORSIKA can only treat primaries with Z=l-26, our cosmic ray spectrum implementation is only 
valid up to ^ 50 PeV. As shown in H, primaries with Z=28-92 as well as an additional ad-hoc component 
are important at energies above ^100 PeV. 

To determine the values of A and 7 corresponding to this model of cosmic rays H, a simulation based on 
CORSIKA version 6.030 was used. The high-energy interaction model QGSJET was shown to be the best to 
describe the muon fluxes observed by AMANDA [7| and several other experiments, and is also the fastest of 
the six high-energy interaction models available in CORSIKA. For these reasons it was used as the default 
model in this work. A detailed comparison with the results of other models is presented in Q. A sample of 
10^° showers with energies above 600 GeV was generated, which took approximately 200 GHz-CPU-days. 

The muon flux formula given above (Equation[l]) was derived assuming a flat atmosphere, and is therefore 
only valid up to zenith angles of less than 70°. In order to describe the full range of zenith angles, 0° — 90°, 
the cos(^^) dependence of the vr and K critical energies valid at zenith angles below 70° can be replaced with 
a 008(6**) dependence as in 0. In this work a convenient parameterization of this substitution is given. 



d^z/dE = ^%E-^^ 
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2 The cos(^) cos{0*) correction 

The muon flux formula as given in the previous section (Equation[T]) cannot be used for 9 > 70°. In [6| the 
critical energies of vr and K, here given simply as E'^^.{9) = 115 GeV/ cos(^) and E^{9) = 850 GeV/ cos(^^), 
were calculated taking into account the curvature of the atmosphere, thus extending the range of applicability 
to all zenith angles. The definition of the critical energy of tt or K, i.e., the energy at which the vr or K decay 
length is equal to the interaction length, is somewhat loose, and depends on the location of the first interaction 
of the cosmic ray primary. Its mean is estimated as =85 g/cm^ in [6|, and the average altitude of muon 
production as ho =17 km (for the direction ^ = 0) in M- In = (1 + 0.04 ■ logio(^/l TeV)) • 85 
g/cm^, i.e., Xjn has only a very weak energy dependence. In this work, values of X„i = Xolog2 =79.6 
g/cm^ and ho =19.28 km are obtained by comparison with CORSIKA-simulated muon flux above 600 GeV 
for the standard US atmosphere parameters (Figure[T]), as shown in Figure|2| 
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Figure 1: Model of the atmosphere. Figure 2: Muon track length. 



Assuming that the interaction length of n and K in air 

pdx ~ p{h{x)) ■ Sxt^^k, 



Top of the atmosphere 



X. 



TT,K 

int 



is approximately constant with energy, and the decay length is proportional 



to energy, l^^, 



, the critical energy (from ^ = 6xTr,K) is 



PTT,^ int ^n,KC 
-C/™ = — 7-r-, — OC 




p{h{x))T^^K PiH^))' 

Therefore, cos(^*) = E^;^ (0) / E^;^ {9) = p{h{x{9)))/ p{h{x{0))). 

At angles 9 < 70° the atmosphere can be approximated as flat and 
x = h/ cos{9). The atmosphere above ^ 200 g/cm^ can be considered isothermal lHJ, therefore the density 
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p{h) oc exp{—h/ho) is proportional to the mass overburden above h, h^piji). The first interaction of a cosmic 
ray primary on average occurs at X 85 g/cm^, where 



X 



p{h{x))dx 



p{h) ■ dh 



h h 



p{h) ■ dh 



hp ■ p{h) 
cos{e) ' 



(2) 



i.e. p{h) oc cos(6'). Therefore in this case cos(6'*) = cos(6'). At larger angles 9* can be approximated as the 
zenith angle of the muon at the point of its production [8|, as indicated in the figure above, since the integral 
in Equation |2l is dominated by the contribution from the atmospheric layer with only a width /iq ~ 6.4 km, 
which can be considered approximately flat. 

The altitude of the first interaction, hi(9), can be determined from X{9) = X(0), where 



X{9) 



Pih)^dh 



p{h) dh 



h,{e) 



sin^(6') 
JT+hJW 



where R is the radius of the Earth. This equation was solved numerically for each 9 and different X = 
X(0). Then cos(6'*) = p{hi{9))/ p{hi{0)) as well as the muon track length from the production point 1(9) = 
^JR^ cos^ (9) + 2Rhi + hf - Rcos{9) can be calculated. The quantities hi{9; X), cos{9*; X), and l{9; X) 
evaluated in this way are then averaged over muon production depth X with weight exp(— X/Xo)/Xo. 
These averaged quantities can be fitted with 
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5-parameter fit to the model (av. dev. 0.04%) 
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Figure 3: Average muon production altitude. 
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Figure 4: cos(6'*) parameterization. 
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model with X„=l 14.8 g/cm', h,|=19.28 m 
3-parameter fit to the model (av. dev. 0.46%) 
5-parameter fit to the model (av. dev. 0.40%) 
model for the South Pole atmosphere 
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Figure 5: Average muon track length. 



Figure 6: Mass overburden of the atmosphere. 



and l{9) 



10^ km 



where x = cos{9). Additionally, the total atmospheric mass overburden as seen under zenith angle 9 was 
parametrized as 

Xtot{9) 



1 mwe 



Pi + P2XP^ + P4(l — x^)P^ 

All fits were first performed assuming the last two parameters are zero, and then another time allowing all 
five parameters to vary, as indicated in Figures 130 The five-parameter fits are accurate to within better than 
0.40% and their parameters are summarized in Table [T] 



fit 


Pi 


P2 


P3 


Pa 


Pb 


av. deviation 


h^i9) 


0.00851164 


0.124534 


0.059761 


2.32876 


19.279 


0.04% 


cos{9*) 


0.102573 


-0.068287 


0.958633 


0.0407253 


0.817285 


0.20% 


m 


1.3144 


50.2813 


1.33545 


0.252313 


41.0344 


0.40% 




-0.017326 


0.114236 


1.15043 


0.0200854 


1.16714 


0.20% 



Table 1: Parameters of the fits to the quantities hi{9), cos{9*), l{9), and Xtot{9) of Section|2l 



The quantity cos(6'*) was also calculated as the zenith angle of the muon direction at the point of its 
production. As shown in Figure |7J the result is almost identical to the density ratio consideration. Also 
shown for comparison are the parameterizations of |5 1. 

3 Muon energy losses and decay 

The muon integral flux at the muon production point i is 
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Figure 7: cos (6'*) parameterization. 



Figure 8: Muon energy loss parameters. 



where dNi/dE is described by Equation [T] with the cos(6') correction of Section|2l The muon integral flux at 
the surface, ff(Ef), can be calculated from 

oo oo 

fAEf) = I dfj{E) = I {dU{El{E))) = (MElm)), 

where E^{Ef) is determined by a stochastic process of muon propagation and depends on a random variable 
C,. To determine the function ff{Ef) one must find the average of fi{E^{Ef)) with respect to all possible 
muon initial energies E^ which result in a muon final energy of Ej. 

On average, final energy of Ej corresponds to some {E^) = Ej + AE. Using the usual linear approxima- 
tion of the muon energy losses dE/dX = a + bE, 

{Ef) = ((a + bEf) exp(6X) - a)/b, (3) 

where X = Xtot{6) — Xq is the average mass overburden the muon crosses before reaching the surface. 
Assuming that the muon energy losses are small, the deviation of E^ from its average {E^) can also be 
estimated. Using an approximation {{E^ — {E^))"^) = (a* + b*E)'^dX on a small path dX, and assuming that 
such deviations can be added on a path X, one derives 

{{Ef - {Ef)r) ^l{a* + b^EfdX = I ^^-^^dE. (4) 

Xf Ef 

Formulae El and m were used to fit a Monte Carlo simulation of muons propagating with energies ~ 600 GeV 
— 60 TeV through 10 and 360 mwe of air performed with MMC [9|. Results of the fits are shown in Figure 
[51 Dashed blue lines show the result of muon propagation through 10 mwe, and solid green lines through 
360 mwe, both rescaled to 1 mwe. The region of the fit is shown with vertical lines and the fits themselves 
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with solid light blue lines. A single set of parameters a* and h* seems to describe well the — {eI)Y) 
deviations for both 10 and 360 mwe thus justifying the assumption of their additivity. 
It is now possible to derive the expression for ff{Ef): 

ff{Ef) = (MEfiEf))) = {M{El) + {El - {Elm 
^ MiEf)) + l:^{{El)){{El - {Elm 



2dE^ 



Therefore, to account for the muon decay it is sufficient to replace Ei in the formula from Section [T] with 
Ej{Ef). Additionally, dEi {dEi / dEf)dEf. All quantities of this section (including the derivative 
dEj/dEf) can easily be evaluated and are implemented in the function fitted to the muon spectrum. 
Muon decay is evaluated as a correction factor exp(— /yuc/ri?/). 

4 Fits to the muon energy spectrum at the surface 

The fit to Formula [T] is shown in Figure 13 and to the corrected formula in Figure [TOl Both figures show 
parameters of the fit A and 7 as well as the of the fit in the region of cos(6') = 0.3 — 1 and in the full 
range cos (6') = — 1. Figure [TT] shows the behavior of the parameters and of the fit as a function of 
the lower boundary of the fit, the upper boundary always at 008(6*) = 1. Equation [T] can only be used for 
cos(6') > 0.3, as mentioned before, and gives the same result in this region as the cos(6')-only corrected 
formula. Both approaches fail to produce values of the fit parameters independent of the region of the fit. 
Including muon energy losses substantially improves the fit. Additionally correcting for muon decay gives 
the best description of the muon spectra. The value of drops to 1.2 at lower boundary of cos{9i) = 0.03 
and is 1.12 at cos(^;) = 0.3. The values of A and 7 vary very little as the region of the fit changes, and are 
0.7015 and 2.715 at cosC^,) = 0.3. 
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Figure 9: Fit with Formula [B to the CORSIKA- 
simulated muon energy spectrum. 



Figure 10: Fit with the corrected Formula [T] 
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The resulting muon differential flux at 10 TeV is plotted in Figure [121 which can be compared with a 
figure containing results of different approaches to calculation of cos(^^*) from 0. Resulting fits to the muon 
integral fluxes are shown in Figures [O] and [HI CORSIKA-simulated muon energy spectra appear to have 
some unexpected structure at cos(^^) < 0.03, which causes a high ^ 50 of the fit to the whole region of 
cos(6') = — 1. This was noticed before fW\, but was shown not to affect the simulation of the muon flux at 
the depth of AMANDA-II (1730 m) or deeper. The highest value of the zenith angle at the surface which can 
be seen at the depth of AMANDA-II is 88.7°, i.e. cos(^^) > 0.023, where the artifact is much less pronounced. 
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Figure 1 1 : Behavior of the parameters of the fit. 



Figure 12: Muon differential flux at 10 TeV. 



5 Uncertainties in the fit model 

To investigate the effect of the uncertainties in the knowledge of the atmospheric profile, the standard US 
atmosphere used in the fits of Section[21was replaced with the winter (July) atmosphere at the South Pole (see 
Figure[l]for comparison with standard US atmosphere). The ground level for this South Pole atmosphere was 
chosen at 1 15.5 m below the see level to match the atmospheric pressure at ground with that of the Standard 
US atmosphere. As seen from Figures [SE only the average muon production altitude and the muon track 
length undergo substantial changes. As seen from Figure [H the change in the average muon production does 
not affect the cos (6**) calculation. Most of the muon decays occur for muons with the longest tracks, which 
come from near the horizon. There the deviation in the track length from that calculated for the standard US 
atmosphere is small, as seen from Figure[51 Therefore, quantities parameterized in Section[21for the standard 
US atmosphere can be used for the atmosphere with quite different profile, in the extreme case the winter 
atmosphere at the South Pole. 

Fits of the previous section were calculated to the CORSIKA-simulated muon data generated in the ab- 
sence of the magnetic field. As seen from Figure [T^ adding a 50 nT magnetic field changes the muon 
scattering profile. The scattering width is 2-4 times that without the magnetic field, but is still very small for 
muons with energies above 600 GeV. 

Fit parameters for both a wrong atmospheric profile and a non-zero magnetic field are shown in FigurelTSl 
As expected, the change from the calculation of Section [4j is negligible. However, the parameters do change 
somewhat when the fits are recalculated in the energy range 1-100 TeV. 
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Figure 13: Integrated muon flux at 600 GeV. 



Figure 14: Integrated muon flux at 10 TeV. 
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Figure 15: Uncertainties in the fit parameters. 



Figure 16: Magnetic field influence. 



6 Energy spectra of atmospheric muon and electron neutrinos 

Along with muons, CORSIKA generates muon and electron neutrinos, spectra of which were fitted with 
the formulae from ifTTl . Muon neutrino spectra were fitted with 

dN 2.85 ■ 10-2 ^ f E^\-^ [ 1 0.213 



dE cm2 sr s GeV V GeV / I i i 6^.^ cos(e*) lAAE.^cosje*) 

\ 121 GeV ~^ 897 GcV 
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and electron neutrino spectra were fitted with 



dN _ 2.4-10-3 
dE cm^ sr s GeV 



0.05 



897 GeV 



+ 



0.185 



1 + 



cos(e*) 
194 GeV 



+ 



1 



1.21E^^ cos(6'*) 
121 GeV 



where ( = a + b logio 



0.11 - 2.4cos(e) 



-0.22 + 0.69 cos(e) cos(^) < 0.3 
-0.46 - 0.54 cos(^) ' " ~ \ -0.01 + 0.01 cos{9) cos(^) ^ 0.3 ' 

The muon neutrino energy spectrum fit demonstrates the same stability with respect to changes in the 
lower cos(6'i) boundary of the fit as the muon energy spectrum fit. Electron neutrino fit parameters, however, 
change in the range A = 0.732 - 0.998 and 7 = 2.689 - 2.732 as 005(61) is varied from to 0.5. The fits 
shown in Figures [TT] and [TS] were produced with cos(6';) = 0.3. 
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Figure 17: Muon neutrino energy spectrum fit. 



Figure 18: Electron neutrino energy spectrum fit. 



7 Results and Conclusions 

Muon and muon- and electron neutrino energy spectra produced with CORSIKA in the energy range 600 
GeV — 60 TeV were fitted with Formula [T] after a cos(6') correction as described in Section |2l (see Table|2l). 
Muon spectra fits were additionally corrected for the muon energy losses and decay. These corrections were 
shown to be important in the selected energy range. Convenient parameterizations of cos(6') correction and 
muon energy losses were given. 













A 


0.701 ±0.006 0.646 ±0.003 0.828 ±0.155 


0.340/0.367 


0.352/0.310 


0.465/0.472 


7 


2.715 ±0.001 2.684 ±0.001 2.710 ±0.023 


2.720/2.712 


2.680/2.696 


2.719/2.741 


A 


with 7 as given above 


0.330/0.374 


0.362/0.285 


0.437/0.382 



Table 2: Muon and muon- and electron neutrino energy spectrum parameters. 
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